////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  ------------------------------ Ebrahim Foulaadvand, 13 May 2012 ------------------------------------- //
//                                                                                                        //
//  The programme "PopulationMaster" evaluates the temporal evolution of the average population number    //   
//  which is subjected to linear birth and death processes.                                               //                                                                                                        //
//                                                                                                        //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>


using namespace std;


double  tau=0.01, beta=0.2,gamma=0.1,Avnum; // tau=time step, beta=birth rate, gamma=death rate.

int M=100,m=20,T=1000,t;

// M=the population number at current time step, m=initial number of individuals, T=number of time steps.

main()

{

ofstream  fileAvnum ("Avn-t beta=0.2 gamma=0.1.plt"); //output file for computed time dependence of pupulation number.
ofstream  fileAvnumanal ("Avnanal-t beta=0.2 gamma=0.1.plt"); //output file for analytic time dependence of pupulation number.

vector<double> P(M+2,0),Avn(T+1,0);

// Array "P" stores the probabilites of having a population number at the current time step.
// Array "Avn" stores the average polulation number at each time step.


// ---------------------------------- initial conditions --------------------------------------------

for(int n=0;n<=M;n++){ P[n]=0; }

P[m]=1;

Avn[0]=m;

//---------------------------------------- Euler Method  -----------------------------------------------

for(int k=1;k<=T;k++){

Avnum=0;

P[0]=P[0] + tau*gamma*P[1]; //cout<<"P[0]= "<<P[0]<<endl;//getch();


for(int n=1;n<M;n++){

P[n]=P[n] + tau*(beta*(n-1)*P[n-1] + gamma*(n+1)*P[n+1] - (beta+gamma)*n*P[n]); //cout<<"P["<<n<<"]= "<<P[n]<<endl;//getch();

Avnum+=n*P[n];

}

P[M]=P[M] + tau*(beta*(M-1)*P[M-1] + gamma*(M+1)*P[M] - (beta+gamma)*M*P[M]);

//cout<<"P["<<m<<"]= "<<P[m]<<endl;getch();


Avnum+=M*P[M];


fileAvnum<<k*tau<<" "<<Avnum<<endl;

fileAvnumanal<<k*tau<<" "<<m*exp((beta-gamma)*k*tau)<<endl;


//cout<<"Avnum= "<<Avnum<<endl;//getch();

//cout<<"Avnumanal= "<<m*exp((beta-gamma)*k*tau)<<endl;getch();

} // end k 

cout<<"end."<<endl;getch();

return 0;

}